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ABSTRACT 

Wc have numerically studied the thermal evolution of various-mass terrestrial planets in hab- 
itable zones, focusing on duration of dynamo activity to generate their intrinsic magnetic hclds, 
which may be one of key factors in habitability on the planets. In particular, we are concerned 
with super-Earths, observations of which are rapidly developing. We calculated evolution of 
temperature distributions in planetary interior, using Vinet equations of state, Arrhenius-type 
formula for mantle viscosity, and the astrophysical mixing length theory for convective heat 
transfer modified for mantle convection. After calibrating the model with terrestrial planets in 
the Solar system, we apply it for O.l-IOM^ rocky planets with surface temperature of 300 K 
(in habitable zones) and the Earth-like compositions. With the criterion for heat flux at the 
CMB (core-mantle boundary) , the lifetime of the magnetic fields is evaluated from the calculated 
thermal evolution. We found that the lifetime slowly increases with the planetary mass (M p ) 
independent of initial temperature gap at the core- mantle boundary (ATcmb) but beyond a crit- 
ical value M CjP (~ O(l)M^) it abruptly declines by the mantle viscosity enhancement due to the 
pressure effect. We derived M CjP as a function of ATcmb and a rheological parameter (activa- 
tion volume, V*). Thus, the magnetic field lifetime of super-Earths with M p > M PjC sensitively 
depends on ATcmb, which reflects planetary accretion, and V* , which has uncertainty at very 
high pressure. More advanced high-pressure experiments and first-principle simulation as well as 
planetary accretion simulation are needed to discuss habitability of super-Earths. 

Subject headings: TERRESTRIAL PLANETS, THERMAL EVOLUTION, MAGNETIC FIELD, 

Introduction recently several super-Earths with masses of a few 

to ten Mq have been discovere d by improved ra - 
Many of exoplanets so far detected may be gas dial velocity measurements (e.g. , llJdrv et alJ2007h 



giants with masses > 100M ffi , because massive Qr micro i cnsin g observations (e.g.. lBeaulieu"et1d 

planets are more easily to be detected. However 
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2006 ). On-going radia l velocity dMavor et al 
20091 ) and microlensing (IGould et al. 201C | sur- 



veys and th eoretical studies (e.g.. llda fc Lin 



2004L 



2008 , 2010l) strongly suggest ubiquity of super- 
Earths in extra solar planetary systems. 

Space transit surveys such as CoRoT and 
Kepler will also detect many super-Earths. In 
fact, CoRoT have detected t he minimum mass 
transiting planet ( CoRoT-7b) (jLeger et al. 20091: 



Queloz et al.l 120091) . With an assumed compo- 



sition, a mass-radius relationship of the super- 
Earths gives planetary masses from transit obser- 
vational data. On the other hand, if the planetary 
masses are obtained by follow-up radial velocity 
observations, the mass-radius relationship can be 
used to estimate the planetary composition, al- 
though there is some ambiguity d e pending on the 
amoun t of H2O ( Sotin et al.ll2007l ). Valencia et al 



(l200fil ) used Birch-Murnaghan equation of state 
(EOS) for rocks and metals to obtain a mass- 
radius relationship of various-mass terrestrial 
planets under some condition s (core rati o , sur - 



face temperature, and etc.). ISotin et al.1 (|2007l ) 
also considered ocean planets which contain 50% 
of H2O. They used the EOS including thermal 
pressure to describe P — V — T relations for ices 
under extremely high pressure and obtained a 
mass-radius relationship for both terrestrial and 
ocean planets. 

Because super-Earths should exist also in hab- 
itable zones, the aspects related to the habitabil- 
ity of super-Earths are being discussed. Plane- 
tary habitability is often discussed in terms of the 
stability of liquid wa ter on the planetary surface 
( Kasting et al.lll993 ). Assuming planets that are 
massive enough to maintain dense atmosphere, a 
range of orbital radius in which liquid water is sta- 
ble is called a "habitable zone." 

In addition to the existence of liquid water, evo- 
lution of amount and composition of planetary 
atmosphere may also be an important factor for 
habitability. It is believed that most fraction of 
the present atmosphere of t he Earth was formed 
by impact degassing (e.g., lAbe fc Matsul Il985l ) 
and it consisted of CO2 and H2O with more than 
100 bars. The plate tectonics on the Earth has 
removed huge amount of CO2 from the Earth's 
atmo sphere on Gyr timcscales (jTaiika fc Matsui 
19921 ). 



(|2007l ) investigated the possibility of plate tec- 
tonics on the surface of super-Earths. The plate 
tectonics would significantly affect amount and 
composition of planetary atmosphere through 
carbonate-silicate cycle with degassing and weath- 
ering. It also has stabilizing effect of planetary sur- 
face temperature, since temperature dependence 
of weathering rate of carbonate works as a nega- 
tive feed back mechanism for su rfa ce temperature 



chang e (jTaiika fc Matsuil 119921) . IValencia et al 
(|2007| ) argued that super-Earths could invoke 
plate tectonics, because terrestrial planets larger 
than Earth may have a thinned surface thermal 
boundary layer and increased yield stre ss. On 
the other hand. lO'Neill fc Lenardk] (|2007f ) showed 
that super-Earths would have stagnant-rid style 
mantle convection without plate tectonics. It is 
noted that these studies a priori assumed ther- 
mal structure of super-Earths, as the studies on a 
mass-radius relation did. 

Although the thermal effect on the mass-radius 
relation may be negligible, plate tectonics should 
depend on planetary thermal evolution history. 
Dynamo activity to generate a magnetic field 
also depends on the thermal evolution. Plane- 
tary magnetic field prevents stellar winds from 
splitting the planetary atmospheres. The mag- 
netic field also prevents cosmic rays from pene- 
trating to the planetary surface. Thus, it may 
be one of the most important factors for land- 
based life to be maintained. It is widely accepted 
that the Earth's ma gnetic field is attributed to a 
dynamo effect (e.g.. iGlatzmaier fc Roberts 19951 
Kuang fc Bloxhamlll997tlKagevama fc Satolll997l ) 



Valencia et al.1 (|2007j ) and lO'Neill fc Lenardk 



in the metallic core. The intrinsic magnetic field 
may be sustained if convective fluid motions in 
the core is vigorous enough, in other words, the 
heat flux through the core surface is large enough. 
Thus, detailed study on thermal evolution of the 
core is needed to evaluate generation of the mag- 
netic field. 

Usi ng the box model (see below) . IPapuc fc Davies 
(j2008l ) modeled the thermal evolution of the 
various-mass terrestrial planets on geological 
timescales to discuss the evolution of planetary 
surface activity, i.e., the plate velocity and the de- 
gassing rate. Their results showed that the super- 
Earths may have dense atmosphere in their early 
history, since larger planets have higher degassing 
rates. Since they were concerned with planetary 



2 



surface activity, they focused on treatment of heat 
generation of radio activity in the mantle and sur- 
face heat flow, leaving treatment of cores simple. 
We will show that careful treatment of the core 
such as the effects of inner core nucleation and 
increase in heat capacity due to high compression 
and gravitationa l energ y stored in the core, which 
Papuc fe Daviesl ( 2008 ) neglected, are important 
for the study of dynamo activity (Note that the 
evaluation of the surface activity is hardly affected 
by the careful treatment of the core). 

In a series of papers, Schubert, Stevenson and 
their colleagues developed a " box" model for ther- 
mal evolutio n of the terrestrial planets in our Solar 
syste m (e.g., Schubert et al.fl~979t Stevenson et al 



19831 ). In their model, the thermal structure is de- 
scribed by two boxes that correspond to the man- 
tle and the core. The temperature variation in 
each box is neglected and the temperature distri- 
butions in the Earth are represented by three dis- 
tinct temperatures of the core, the mantle, and 
the planetary surface. The heat flow is evalu- 
ated by thermal conduction through the thermal 
boundary layers at CMB and the planetary surface 
with the thermal Boun dary Layer Theory (BLT) 
(jStevenson et al.l Il983j) . In the BLT, the heat 
flux is determined by the thickness of the thermal 
boundary layer, which is given by local Rayleigh 
number. Because this model is easily treated, it 
provides a powerful tool to explore general trends 
of thermal evolution of terrestrial planets. 

For Earth's thermal evolution, we can use ob- 
servational constraints such as surface heat flow 
and inner core size. With the calibrated model, 
unknown parameters such as ini tial temperature 
distribution of Earth's interior ( Yukutake 2000l 
also see section 2.7) or pota ssium abundance in 
the core ( Nimmo et al. I l2004h can be constrained. 
Rhcological parameters and impurity abundance 
in the core are also estimated (see section 3.3). 
The existence of the magnetic field for Mercury 
and early decay of the magnetic field for Venus and 
Mars are consistent with calculations with rea- 
sonable choice of initial tempera ture or impurity 



abun dance in the cores by BLT (jStevenson et al 
1983h and MLT (Appendix B) 



Gaidos et al.l ( 2010l ) simulated thermal evolu- 



2.5M®) can not sustain magnetic field, because 
they found an inversion of gradients of the melt- 
ing and adiabatic curves in the core under such 
high pressure. Since the outer core is solid in the 
inverse state, the cooling of the inner liquid core, 
in which dynamo operates, is inhibited. Although 
the possibility of the inversion raised by the paper 
is a very important, the conclusion depends on 
high pressure material properties such as meting 
and adiabatic curves that need to be confirmed. In 
the present paper, we point out another important 
factor to inhibit dynamo activity in super-Earths, 
drastic increase in the mantle viscosity due to pres- 
sure effect. Even if the inversion in the core does 
not occur, the enhanced mantle viscosity quickly 
terminates dynamo activity in the core. 

In the present paper, we are concerned with 
lifetime of magnetic fields of super-Earths. There 
is no observational constraint for super-Earths 
to calibrate parameters for a model and ini- 
tial/boundary conditions. Here, clarifying key 
quantities for generation of the magnetic field, 
we derive planetary mass dependence of lifetime 
of magnetic activity and how it depends on the 
model parameters and initial/boundary condi- 
tions. To reduce unknown parameters, we con- 
sider super-Earths in nearly circular orbits in 
habitable zones. The analysis on the key rhco- 
logical parameters will provide new motivations 
for high pressure experiments and first principle 
simulations, since state-of-arts high pressure ex- 
periments have already reached the pressure at 
the bottom of Earth's mantle. We will point out 
that initial temperature distribution sensitively 
affects the lifetime, which gives new motivations 
for theories of planet accretion from planetesimals 
and core formation. 

Here, we develop a thermal evolution model to 
discuss the existence of the intrinsic magnetic field 
in terrestrial planets with various masses. Since 
we are concerned with heat flux across the CMB, 
we calculate detailed radial temperature distribu- 
tion in both core and mantle. We use the Mixing 
Length Theory (MLT) to calculate the heat flow. 
The MLT is commonly used to study stellar inte- 
rior. We use the modified version for low Reynolds 
number flow in solid planets that have radial dis 



tion of various sized super-Earths by the BLT to 
evaluate the magnetic activity of super-Earths. 
They concluded that massive rocky planets (> 



continuities in their interior ( Sasaki fc Nakazawa 



198fit lAbl[l995l: ISenshu et~ail 



2002t iKimura et al 



20091 and references therein). The modified 
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MLT is useful for calculations of super-Earths that 
may have additional higher-pressure phase tran- 
sitions such as a post-post perovskite transition 
(although we do not consider it in the present 
paper) and early-stage planets that may have 
convectio n-barriers at upper /lower mantle bound- 
ary (e.g.. iHonda et ah 19931) or density cros sover 
at melt / solid boundary Labrosse et al.l ( 20071 ) . We 
will also show the results for the two-layer con- 
vection case in which convective flow does not 
penetrate the spinel-perovskite transition at up- 
per/lower mantle boundary, while most of calcu- 
lations are in the cases of one- layer convection. We 
compare the MLT with the conventional BLT in 
details and show that they are in good agreement 
with each other in the case of one-layer convection 
(Appendix A). 

We will show that the lifetime is rather shorter 
for super-Earths than for Earth-mass planets for 
nominal parameters of solid state. The mecha- 
nism to suppress dynamo activity in super-Earths 
found in thi s pap er is independent of that in 
Gaidos et all (|2010h . so that it may be likely that 
super-Earths are magnetically inactive. We also 
point out that the choice of initial conditions and 
rhcological parameters highly affect the thermal 
evolution of the planets. In section 2, we explain 
our numerical model. The numerical results will 
be shown in section 3. Finally we discuss the hab- 
itability for terrestrial planets in view of the in- 
trinsic magnetic field. 

2. Numerical Model 

We follow thermal evolution of a planet that 
is a cooling process from a hot early state due 
to accretion from planetesimals and core-mantle 
differentiation, for various-mass terrestrial planets 
by using a one-dimensional spherically symmetric 
model. 

As described below, there are unknown param- 
eters for rhcological properties and initial condi- 
tions of planets. Furthermore, even in our So- 
lar system, the terrestrial planets have different 
compositions (a water-rock-iron ratio) and surface 
temperature that is regulated by orbital radius. 
In extrasolar planetary systems, more variations 
in the water-rock-iron ratio and surface temper- 
ature should exist due to different metal abun- 
dance of host molecular clouds as well as differ- 



ent pressure/temperature state of disks and plan- 
etary atmosphere and planetary formation pro- 
cesses. These unknown variations could make the 
study on the mass dependence meaningless. 

In order to reduce the uncertainty, in a "nom- 
inal" case (see section 2.8), we consider planets 
with surface temperature (T sur f) of 300K and the 
mantle/core mass ratio (Cm/c) of 7 : 3, which is the 
same as that for the Earth. The planets may cor- 
respond to extrasolar terrestrial planets in nearly 
circular orbits in habitable zones. In the nomi- 
nal case, the other unknown parameters in mantle 
rheology, core impurity and initial temperature are 
calibrated with data of the Earth. 

The restriction to the nominal case enables us 
to derive clear planetary mass dependences. We 
also discuss how these dependences change by dif- 
ferent choice of the parameters. In Appendix B, 
we also carry out the runs with different T sur f and 
Cm/c t° calculate thermal evolution of Mercury, 
Venus and Mars. The results show that our model 
can be applied for these planets too, if we use rea- 
sonable non-nominal parameters. Systematic sur- 
vey for thermal evolution of super-Earths in the 
non-nominal cases is left for future works. 

The thermal evolution is calculated by the fol- 
lowing methods: 

1. The radial density distribution is calculated 
by using VINET equation of state taking 
into account pressure dependence (section 
2.1). Since this distribution is almost in- 
dependent of evolution of temperature dis- 
tribution and the inner core growth as ex- 
plained below, we use the distribution calcu- 
lated in the initial state (t = 0) throughout 
the entire thermal evolution. 

2. Given a temperature distribution in the in- 
terior at time t, the radius of the inner solid 
core is calculated by the melting tempera- 
ture with pressure and composition depen- 
dences (section 2.2). Sulfur is considered 
as the impurity in the core and its con- 
centration in the outer liquid core is self- 
consistcntly calculated with the condensa- 
tion of the inner core (section 2.2). 

3. The heat transfer throughout the mantle 
is calculated by the astrophysical mixing 
length theory modified for solid planets (we 
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discuss its validity and usefulness in section 
2.3). The mantle viscosity in the heat trans- 
fer equation is estimated by using Arrhcnius 
type formulation (section 2.5). By subtract- 
ing the energy loss during the timestep (At), 
we obtain a new temperature distribution at 
t + At and go back to step 2. Time evolution 
is calculated by iterations of step 2 and 3. 

In the following subsections, we will explain each 
step in more detail. 

2.1. Density profile 

The hydrostatic stratification is calculated by 

dP dr 1 

~dr = P{r)9{r); d^ = (1) 

where P, r, p, g, and m are pressure, radius, den- 
sity, gravitational acceleration, and mass inside 
the radius r, respectively. Vinet EOS (|Vinet et al 



19871 ) is given by 



P = 3A - 



1 



■exp(0(l - x)), 



where 



= 2^0-1). 



(2) 



(3) 



x = po/p, and Po,Kq and K' are density, bulk 
modulus, and its pressure differentiation at zero 
pressure, respectively. 

Birch-Murnaghan EOS is "finite strain" EOS, 
in which pressure is expressed by Taylor series ex- 
pansions of finite strain, and it is often used for the 
calculation of interior of solid planets. However, 
the finite strain EOS do not accurately represents 
the volume variation under very high compres- 
sion (if pressure exceeds the bulk modulus of zero 
pressure, the expansion never converges). Vinet 
EOS is derived from a general inter-atomic po- 
tential energy function. For simple solids, Vinet 
EOS provides more accurate representations of the 
volume variations with pressure, under very high 
pressure. Since we are concerned with interior of 
super-Earths under very high pressure, we adopt 
Vinet EO S rather than Birch -Murnaghan EOS, 
following ((Valencia et al.ll2007l) . 

Since thermal contraction is small enough (it 
is less than 1% in physical size for 100K change 
in the average temperature of Earth's interior), 



we neglect the temperature dependence in Vinet 
EOS. The density profile in the planet is calculated 
by numerically solving equations ([1]) and ([2"]). 

The compositions are assumed as olivine and 
7-spinel for upper mantle, perovskite and post- 
perovskite for lower mantle, Fe and FeS for outer 
core, and Fe for inner core with the properties 
given in Table Q] As temperature decreases, the 
inner solid core grows and sulfur moves from the 
inner core to the outer core (see section |2~2"]) . In 
general, the inner core growth changes the vol- 
ume of the whole core, because the parameters, 
Po, Kg and Kq, are different between Fe and FeS. 
However, since we numerically found that the vol- 
ume change of the whole core is very small, we 
neglect it. Figure [1] shows the result of numeri- 
cal calculation of the density profile for 0.1 to 10 
Earth-mass planets. This result is different in the 
radii of the cor e and t he m antle from the resu lts by 
Valencia et all (|2006l) and lSotin et~aTl ([20071) by a 
few %, which may be due to different choice of 
parameter values in the EOS. But, this difference 
docs not affect the thermal evolution. 

2.2. Thermal evolution of the core 

The temperature distribution of the core is de- 
termined as follows: 

1. The inner solid core: We assume that each 
part of the inner core memorizes the temper- 
ature at which it solidified, because of ineffi- 
cient heat transfer due to conduction in the 
solid core. 

2. The outer liquid core: We assume that the 
liquid core has an adiabatic temperature dis- 
tribution by vigorous convection. 

3. Time evolution: The radius of the inner core 
and the temperature at the CMB is deter- 
mined by total energy of the inner and outer 
cores as described below, and the total en- 
ergy is given as a function of time with inte- 
grating heat flux at the CMB. 

The adiabatic temperat ure gradient i n the 



outer core is given by (ISohl fe Spohn 
Yukutakdl2000l IValencia et al.ll2006l) 



1997 



dT 
dr 



P91G 



(4) 
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where 70 and K s are Griineisen parameter and 
bulk modulus of the liquid core. Depth varia- 
tion of 7g is calculated as = 7Go(po/p) 9 (the 
parameter values used are summarized in Table 
[T|). The density at pressure pooc, bulk modu- 
lus KqoCi and its pressure derivation K' GOC of the 
outer core are given by impurity concentration xs 
as: 

Zfc + Zs 

XFeS = XS - 



Pooc 



1 - ZFeS 
PFc 



PFcS 



K< 



1 



00c 



Pooc 



-xfcS 1 



pFe 



K Fc 



PFeS K?b\ 



K' 



00c 



oc 



1 - X Fe S 1 



K' Fo 



(5) 



(G) 



(7) 



a^FeS 1 



where M- lc and M c are the inner core mass and 
total mass of the inner and outer cores and xqs is 
the initial impurity concentration. In the nominal 
case, we adopt xos =0.1. 

Given the inner core radius, we can calculate 
the total energy of the core (E corc ) which is sum 
of the gravitational energy (E g ), latent heat (£"/), 
and thermal energy (E t h)- As described above, the 
temperature at the CMB is given as a function of 
the radius of the inner core. As a result, we can 
obtain -E^orc as a function of the temperature at 
the CMB. Conversely, the radius of the inner core 
and the temperature at the CMB are given as a 
function of E corc . 

The energies are given by 



PFc 



K 2 

"Fc 



PFcS 

(8) 

where xfc, ^FeS> Zp e , and Zg are mass fraction of 
Fe and FeS, molar weights of Fe and S, respec- 
tively. 

The inner core nucleation decelerates cooling 
of the core by release of gravitational energy 
due to the change in the den sity distribution 



and by release of latent heat (jStevenson et al 



19831 : iGubbins et all 120041 ) . The light elements 
are kicked into the outer core, resulting in de- 
pression of a melting point of the ou ter core 
( Stevenson et al. 1983 ; lYukutakd l2000h . The 
boundary between inner and outer cores is lo- 
cated at the intersection between adiabatic and 
melting curves in the core. We use a Lindeman's 
equation for the melting curve of pure iron, 



t(p) = r n 



2/:! 



exp 



270 

q 



1 - 



(9) 

We also consider the depression of a melting point 
by concentration of light elements. We define the 
melting point of Fe-FeS alloy as 



T me it = (1 - 2x s )T(p), 



(10) 



and the factor (1 — 2x§) expresses the depression 
of the m elting point due to dissolution of ligh t 
elements ( Usselmanlll975t IStevenson et al. 19831) . 
Assuming that the outer core is well mixed by con- 
vection, 

M r . 



xs = xos 



M r - Mu 



(11) 



K 2 

"FcS 



Ei 
Eth 



LM„ 



Airr 2 p{r)C p {r)T{r)dr, 



4iTr 3 p oc {r)g oc (r)dr, 



(12) 



where L is the latent heat released by solidifica- 
tion of unit mass of iron, w hich is assumed to 
be co nstant of 1.2 x 10 6 J/kg (lAnderson fe Duba 



19971 ). and not to depend on the impurity concen- 



tration in the outer core, and C p is specific heat 
with constant pressure. Both the gravitational en- 
ergy and the latent heat are released after the in- 
ner core starts to solidify. Gravitational energy 
is also released by thermal contraction, which will 
be discussed in section 2.6. The total energy E COTe 
decreases with the rate that is equal to the heat 
flux at the bottom of the mantle (see section 2.3). 
Detailed calculations of the energies are given in 
Appendix C. 

2.3. Heat transfer throughout mantle 

The mantle is cooled by irradiation from the 
planetary surface and heated by heat flow from the 
core and internal radioactivity (see below). The 
heat transfer equation is: 



pC p 



dT 
~dt 



1 d 



r 2 k c 



dT 
dr 



r 2 k v 



dT 
dr 



dT 
Or 



-pQ, 



(13) 

where k c is thermal diffusion coefficient, Q is ra- 
dioactive heat production rate, (dT/dr) s is the 
adiabatic temperature gradient, and the first and 
second terms in the right hand side represent con- 
ductive and convective fluxes. 



G 



To evaluate the convective flux in the man- 
tle, we use the astrophysical mix ing; length theory 
(MLT) modified for solid planets ( Sasaki fc Nakazawal 



19861: lAb j[l995t ISenshu et al.ll2002t iKimura et al 



20091 and references therein) , rather than the con- 
vent ional parameterized co nvection model (PCM; 
e.g. JSharpe fc Peltierlll979h or th e commonly used 
boun dary layer theory (BLT; e.g., Stevenson et al.l 
1983). 

The PCM is very simple (Appendix A). How- 
ever, it uses the values of k c and Rayleigh num- 
ber Ra that represents the whole mantle, which 
are difficult to evaluate for real mantle because of 
huge spatial variation of the mantle viscosity. As a 
result, although the PCM can be applied to study 
overall trend of thermal evolution, it may not be 
accurate enough for evaluation of heat flux across 
the CMB (-Fcmb): which we are concerned with 
in the present paper. In the BLT, since the heat 
flux is expressed by quantities only in the thermal 
boundary layer (Appendix A), the BLT has bet- 
ter resolution for evaluation of i*cMB • Since the 
modified MLT also uses local values of physical 
quantities, it quantitatively agrees with the BLT 
for wide range of parameters, while the PCM does 
not agree with the BLT and the MLT for the cases 
in which viscosity variation is large in the man- 
tle, as shown in Appendix A. As explained below, 
since the MLT is more easily to be applied for 
super-Earths, we use the MLT. 

In the early Earth, the upper/lower mantle 
boundar y could have worke d as a barrier for con- 



vection (jHonda et al.lll9930 . Tentative stagnancy 



at the upper/lower mantle boundary is also sug- 
gested for some subduction slab s in the present 
Earth (jWortel fc Spakmanl [2000h . The density 
overturn at melt/solid boundary in deep magma 
ocean in the e arly Earth may have also worked 
as the barrier ( Labrosse et al ] l2007h . In super- 
Earths, post-post perovskite transition in deep 
man tle at high pressure co uld also work as a bar- 



rier (jUmemoto et al.l l2006) 



As explained below, the modified MLT is eas- 
ily applied for mantle convection with barriers, 
without tuning of parameters for each barrier. Al- 
though most of our calculations in the present pa- 
per only consider the surface boundary and CMB 
(in some runs we consider the upper/lower mantle 
boundary at spincl-pcrovskitc transition as well), 
wc use the modified MLT for future extensions of 



calculations with various convection barriers. In 
the following papers, we will consider the effects 
of other barriers. 

In the MLT, the coefficient for convective heat 
transfer is given by 



p 2 C p ag£ A 





dT 
dr 



for 



dT 
dr 



dT 
dr 



dr 



where rj and £ are viscosity and the mixing length, 
respectively. Here, the velocity of fluid blobs is 
evaluated by Stokes velocity rather than free fall 
velocity in the original MLT, in order to apply 
the model to low Reynolds number flow in the 
mantle. In the astrophysical context such as stel- 
lar interior, the density scale height is usually 
adopted as £. For calculation of thermal evo- 
lution of the Earth, it is proposed that a dis- 
tance (D) from the closest barrier such as the 



priate for £ ( 


Sasaki & Nakazawal 1986; lAbd 19951: 


Senshu et al. 


20021 Kimura et al.l 20091 and ref- 



erences therein). The detailed comparison with 
the PCM and BLT in Appendix A shows that 
£ = 0.82£> is the best choice. We adopt £ = 0.82£> 
for all the runs in the present paper. 

With this choice, as approaching a barrier, k v 
rapidly decreases in proportion to £ 4 and the con- 
ductive term dominates in Eq. (|T3|) . As a result, 
thermal boundary layers, in which the conductive 
heat transfer dominates, are automatically repre- 
sented. Thereby, the modified MLT is easily ap- 
plied for calculation for thermal evolution of the 
proto-Earth or super-Earths. 

The mantle viscosity is discussed in section 2.5. 
The heat flux and temperature at the CMB de- 
termined from the calculation of the mantle heat 
transfer is used as a boundary condition for ther- 
mal evolution of the mantle. The total energy in 
the core is interpolated from the result of section 
2.2. It decreases with time, according to the cal- 
culated heat flux at the CMB. 

2.4. Internal heat source 

For thermal evolution of terrestrial planets on 
geological timescales, long-lived radiogenic ele- 
ments ( 40 K, 232 Th, 235 U, and 238 U) are impor- 
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tant heat sources in the mantle. The estimated 
amounts of these elements in the E arth arc com- 
piled in table [2] ( van Schmusl 1995 ). Here we as- 
sume the same abundances of the radiogenic ele- 
ments in the mantle in super-Earths as those in 
the Earth and that the elements are distributed 
uniformly throughout the mantle. The heat pro- 
duction rate at time t, Q(t), is given by Q(t) = 
HUexp(—X(t — where U, H, £ ffi and A are 
the abundance, heat production rate, and decay 
constant of the clement, respectively. 

2.5. Temperature and pressure depen- 
dency of mantle viscosity 

The mantle viscosity is one of the most im- 
portant physical parameters to simulate thermal 
evolution, since it determines the heat trans- 
fer efficiency in the mantle (sec Eqs. [TJ] and 
P3]). The viscosity sensitively depends on tem- 
perature and pressure, and both of temperature 
and pressure widely vary throughout the man- 
tle. Here, we adopt Arrhcnius type formulation 
for tem perature- and pressure-dependent viscosity 
model ( Ranalli|[200lh : 



ri{T,P) = \ 



1 



■ exp 



E* + PV* 
nRT 



;(!-«)/« 



(15) 

where R, e, n, B, E* and V* are universal gas con- 
stant, strain rate, creep index, Barger coefficient, 
activation energy, and activation volume of man- 
tle, respectively. We use different values of these 
parameters for upper and lower mantles. The min- 
eral properties we use are listed in Table [3] Note 
that the prescription for the mantle viscosity may 
include uncertainty. The formula is based on the 
theoretical rate equation for creep law of rocks. 
In this formula, the most important parameter to 
study thermal evolution of super-Earths is the ac- 
tivation volume (V*), since V* determines the de- 
pendence of the viscosity on pressure and the pres- 
sure in the mantle can be increased by orders of 
magnitude as the planetary mass increases. The 
activation volume is related to atomic volume, but 
the exact values under extremely high pressure is 
not well determined. In the nominal case, we use 
V* = 10 x 10~ 6 m 3 mol _1 , but we also test a smaller 
value of V* = 3 x 10~ 6 m 3 mol _1 . In section 3.4, 
we will discuss how the conclusion in the present 
paper depends on a formula for the viscosity. 



2.6. Release of gravitational energy by 
thermal contraction 

Although the thermal contraction is negligible 
for physical radius, the gravitational energy re- 
leased by the thermal contraction cannot be ne- 
glected (it is about 50kJ/kg for 100K change of 
the core in an Earth- mass planet). In our model, 
the released energy is re garded as increa se in the 
specific heat of the core (lYukutakdl2000h : 



AC P = — . 



aP 
P 



(16) 



The gravitational energy released by thermal con- 
traction is more effective in deeper regions. For an 
Earth-mass planet, AC p /C p is as large as 50% at 
the CMB in our calculation. 

2.7. Initial conditions 

Initial temperature distribution in th e mantle is 
determ ined by the procedure following lYukutake 
(|2000h . as illustrated in Fig. H (1) An adiabat 
is drawn from the bottom of the surface bound- 
ary layer with 1500 K down to the top of the 
boundary layer at the CMB (the obtained ten- 
tative temperature is denoted by T%), assuming 
efficient thermal convection, (2) the initial tem- 
perature at the bottom of the CMB is assumed to 
be Tcmb = T 2 + AT C mb, where AT C mb is deter- 
mined by step 6, (3) the adiabat is drawn from the 
CMB with temperature Tqmb to the surface, (4) 
the initial temperature distribution in the man- 
tle is given by the average of the two adiabats 
obtained by steps 1 and 3, (5) core temperature 
is determined by the procedure given in section 
2.2 with Tomb, and (6) the amount of ATcmb 
(~ 1000K) is determined by the requirement that 
the predicted surface heat flux and inner core ra- 
dius for an Earth-mass planet are comparable to 
the observed values for the present Earth. We use 
this value for all the cases with various planetary 
masses. Note that the temperature distribution is 
quickly relaxed to an equilibrium distribution, as 
long as we use the initial conditions created by the 
above procedures. 

2.8. Simulation parameters 

We summarize parameters for the "nominal" 
case: 
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• Boundary conditions 



3. Numerical Results 



— surface temperature: T sur f = 300K 

— a mass ratio between mantle and core: 
Cm/c = 7:3 

— the CMB is a barrier for convection, 
while convection penetrates the up- 
per/lower mantle boundary 

• Initial conditions 

— impurity fraction: xqs = 0.1 

• Rheological conditions 

— activation volume: V* = 10xl0 _6 m 3 mol _1 

We first investigate planetary mass depen- 
dence for planets with the above nominal pa- 
rameters. For Earth-mass planets, we adopt 
ATcmb = 1000K in most runs, because the 
nominal case with AT CMB = 1000K reproduces 
the present Earth. We also systematically study 
dependences of the results on AXcmb, because 
ATcmb is not well determined. In some runs, the 
upper/lower mantle boundary is treated bar- 
rier for convection. We also carry out calculations 
with different values of T SU rf, Cm/c, and Xqs to re- 
produce the results that are consistent with the 
current magnetic fields of Mercury, Venus, and 
Mars in Appendix B (we do not systematically 
survey the dependences on these parameters). 

2.9. Definition of lifetime of planetary in- 
trinsic magnetic field 

To drive dynamo action, liquid metallic core 
must be in act i ve con vection state. Following 
Stevenson et al.l (1983), we adopt the threshold 
heat flux in the core for generation of dynamo ac- 
tion (conducted heat flux along the core adiabatic 
thermal structure) as, 



-Fcrit = k 



fdT { 



CMB 



\ dr 



= k c P9 J G T C mb- (17) 



We define the lifetime of magnetic field as a pe- 
riod during which the core heat flux exceeds the 
threshold value. 



3.1. Thermal evolution of the Earth 

We now show the evolution of temperature dis- 
tribution calculated by the procedures in section 2. 
Figure^ shows the evolution of thermal structure 
of an Earth-mass (M = 1M©) planet for the nom- 
inal case. Cooling of the mantle slows down with 
time, since the decrease in temperature enhances 
the mantle viscosity (eq. [15] ) and hence depresses 
the efficiency of heat transfer in the mantle. This 
implies that the initial mantle temperature dis- 
tribution hardly affects the thermal evolution on 
timescales longer than Gyr, as long as the ini 



tial t emperature is high enough (jStevenson et al 
1983|). However, since the core works as a heat 



bath for the mantle, the initial Tqmb would affect 
the thermal evolution of the mantle. 

Figures |3Jd, c and d show the time evolution 
of the surface heat flux (-F SU rf), the heat flux 
across the CMB (-Fcmb) and the inner core ra- 
dius (i?i c ), respectively, in the nominal case. We 
adopt ATcmb = 1000K. The inner core emerges 
at 2.2 Gyr. Its radius reaches 1200km at 4.5Gyr 
that agrees with the observed value of the present 
Earth, as expected. With higher/lower values 
of Tcmb, the inner core radius at 4.5Gyr is 
smaller/larger. The growth rate decreases with 
time, because not only geometric effect but also 
the increase in impurity concentration in the outer 
core depresses the melting point of the outer core. 
After the emergence of the inner core, the heat flux 
reduction becomes more slowly because the inner 
core growth releases gravitational energy and la- 
tent heat that work as internal heat sources. The 
heat flux remains larger than the critical value 
given by Eq. (fT7| . which is expressed by the dot- 
dashed line in Fig. [3ji, during first 12 Gyr and the 
lifetime of dynamo activity is expected to be 12 
Gyr for this nominal case. 

Some paleomagnetism data suggests that the 
magnetic field of the Earth i s enhanced to be the 
present level at ~ 2 Gyr (e.g.. lflalelll987l ). It might 
be due to formation of the inner core because nu- 
cleation of the inner core provides additional heat 
source. Paleomagnetism data may include a large 
uncertainty. If more detailed data is provided, it 
will constrain the condition of generation of intrin- 
sic magnetic field. 
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Figures [4] and [5] show the results of the two- 
layer convection. In this calculation, we set the 
upper/lower mantle boundary as a barrier for con- 
vection. The other boundary conditions and the 
model parameters are the same as those in the 
nominal case. The mixing length is shorter in 
the entire regions of mantle and cooling is slower 
than in the one-layer case. As is shown in Fig- 
ures HI if we adopt initial ATcmb = 1000K as 
in the case of one-layer convection, the inner core 
can not grow to 1200km because of the low heat 
transfer efficiency of the layered convection and 
core temperature is somewhat higher than that 
obtained by the one-layer convection. We also car- 
ried out a calculation with initial ATcmb = 800K. 
Figure [5^ shows the evolution of thermal profile. 
A thermal boundary layer at upper/lower man- 
tle boundary is clearly established. Because the 
lower initial ATcmb is compensated with the inef- 
ficient two-layer convection, evolution of heat flux 
through CMB (Tcmb) and lifetime of magnetic 
field (llGyr in this case) are similar to those in 
the one-layer convection case. 

If the two-layer convection is assumed only in 
the Archean and Hadean (t < 2Gyr), with the 
same initial ATcmb (1000K) the surface heat 
flux at 4.5Gyr is T surf ~ 0.12Wm~ 2 , which 
is somewhat higher than the observed value 
(~ 0.09Wm~ 2 ), since the thermal energy beneath 
upper/lower mantle boundary have been stored 
until 2Gyr and then supplied to upper mantle af- 
ter 2Gyr. However, the evolution of Tomb is not 
so different between one- and two-layer convec- 
tion. The lifetime of magnetic field is about 12 
Gyrs. 

Thus, if we tune the initial ATcmb with the 
present observed values of T sur f and i?i c , the ex- 
pected lifetime of magnetic field is not affected by 
the mode (one-layer or two-layer) of mantle con- 
vection. 

3.2. Thermal evolution of Mercury, Venus, 
and Mars 

The existence of the magnetic field for Mercury 
and early decay of the ma gnetic field for Venus an d 

" (Il983h . 



Mars were addressed by IStevenson et al 



using the box model with different parameter val- 
ues such as surface temperature and a mantle-core 
mass ratio from those in the nominal case. The 
validity of these parameter values is discussed in 



Appe ndix B. Adopting the same parameter val- 
ues as IStevenson et al. ( 1983 ) , we have performed 
simulations for Mercury, Venus, and Mars with 
our model. As discussed in Appendix B, our 
model produces the results that are not inconsis- 
tent with the magnetic activity of Mercury, Venus, 
and Mars. 

3.3. Thermal evolution of super-Earths 

For super-Earths, we use the nominal param- 
eters (the surface temperature is 300K and the 
mantle/core mass ratio is 7:3), assuming that their 
orbits are nearly circular and in habitable zones. 
We also assume one-layer convection throughout 
the mantle. Detailed study on the effects of phase 
transitions is left to future works. Figure UK shows 
evolution of temperature distribution for a planet 
with mass M p = 5M®. Compared with the case 
of M p = IMq in Figure a thicker thermal 
boundary layer is established on the CMB within 
first few Gyrs, since the viscosity of the bottom of 
the mantle is higher. The increase in the viscos- 
ity due to higher pressure dominates the decrease 
due to higher temperature (Eq. [IS]). Thus, Tcmb 
is lower than the case of M p = 1M§ (Figs. HJd 
and [3Jd). On the other hand, the effective heat 
capacity rapidly increases with M p . Figures [15] 
in Appendix C show that for fixed Tcmb, ther- 
mal energy T tri oc M 2 , while the core surface area 
S^ore increases with M p only weakly (Score oc M p ). 
Therefore, the core for higher M p cools much more 
slowly. It is shown that R lc docs not grow at all for 
lOGyr. On the other hand, T sur f is not so different 
from that in the case of M p = IMq. The heat bath 
of surface heat flow is radiogenic elements in the 
mantle and that is proportional to M p . To balance 
heat generation and cooling, T surr should be pro- 

1/3 1/3 

portional to M p , provided R p oc M p . Thereby 
T sur f changes by a factor of only 5 1 / 3 ~ 1.7. 

As discussed in the above, thermal evolution 
of super-Earths differs from that of Earth-mass 
planets in many aspects. Here, we focus on evo- 
lution of heat flux through CMB, Tcmb, and in- 
ner core radius, _R; C , in order to study magnetic 
activity of super-Earths. Figures [7] show the evo- 
lution of Tcmb (left column) and R lc (right col- 
umn) with various initial ATcmb for the case of 
(a) M p = lM ffl> (b) 2A/ ffi , (c) 5M e , and (d) 10M e . 
Solid, dot, dashed, and long-dashed lines repre- 
sent the results with initial ATcmb = 1000K, 
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2000K, 5000K and 10000K, respectively. In all 
cases, V* = 10 x 10~ 6 m 3 /mol. 

The results in the left column show that Fcmb 
is generally higher for higher initial ATcmb- The 
dependence is more pronounced for relatively large 
M p cases. For M p = 1 M®, the dependence 
is very weak. We found the dependence is also 
very weak for M p < 1 M®. For M p < 1 M®, 
the temperature dependence of the mantle vis- 
cosity (Eq. [15]) dominates over the pressure de- 
pendence. Then, Fcmb is high when the core 
temperature is high, and -Fcmb declines as the 
core cools. Thus, the heat flux is self-regulated 
to be quickly relaxed independent of the initial 
values. On the other hand, as will be shown 
later, when ATcmb ~ 1 000 (M p /M®), the pres- 
sure dependence is more effective. Then, the self- 
regulation does not work and the dependence of 
-Fcmb on initial ATcmb is retained for more than 
20 Gyrs. The threshold flux for driving dynamo 
action is marked by an black lines in each case. 
The decline of the threshold value is due to de- 
crease of core surface temperature (Eq. [E]). The 
duration for Fcmb > -Fcrit determines lifetime of 
magnetic field generation. 

Papuc k, Daviesl (|2008l ) obtained Fcmb °c 



2 /3 

M p ' , whereas our results shows Fcmb oc 
provided that ATcmb is sufficiently high, 
difference may come fro m the assignment of spe- 



Mp 
The 



cific heat of the core. Papuc fc Daviesl ( 20081 ) 



assumed constant specific heat of the core, C p = 
lOOOJkgKT 1 for all sized planets. As we discussed 
in section 2.6, however, thermal contraction re- 
sults in increase in the effective C v and the effect 
is more pronounced for larger M p . In our calcula- 
tions that include this effect, the core tends to cool 
less efficiently and the dependenc y of Fcmb on M„ 



i s stro nger than that obtained bv lPapuc fc Daviea 
(120081 ). 

The right column shows the growth of inner 
solid cores for AT CMB = 1000,2000,5000 and 
10000K. In the case of M p = 1 M®, an inner 
core is nucleated at 2-3 Gyrs, almost indepen- 
dent of initial ATcmb, since the core cooling is 
self-regulated. For M p = 2M®, the inner core 
growth depends on ATcmb for ATcmb ~ 2000K. 
For such high ATcmb, since the core has larger 
thermal energy initially and the heat flux is not 
self-regulated, it takes more time for the core tem- 
perature to become below the nucleation temper- 



ature. For M p > 5M®, core hardly cools on 20 
Gyrs, the inner core does not grow from the ini- 
tial state. In these cases the inner core size is 
determined by a relationship between adiabatic 
curve and melting curve of iron. The inner core 
of super-Earths (M p > M®) have never nucleated 
for ATcmb = 10000K. For massive planets, the 
increase in the viscosity due to higher pressure is 
overcome only by very high initial temperature. 
The high ATcmb also delays nucleation of the in- 
ner solid core. As a result, there is trade off be- 
tween heat flux and inner core growth through the 
relation between melting point of iron core and 
temperature- and pressure-dependency of mantle 
viscosity. 

Figures [8] show dependence of evolution of 
Fomb on M p for fixed values of ATcmb- For 
ATcmb = 1000K, we have already mentioned 
that Fcmb is rather lower for M p = 5M® than 
for M p = M® (Figs. |Bb and [5b), because the 
increase in the viscosity due to higher pressure 
dominates the decrease due to higher temperature 
for M p = 5M®. This trend is clearly shown in 
Fig. Et. 

However, this is not always the case. If the 
core temperature is high enough (in other words, 
ATcmb is high enough), or if pressure is low 
enough (M p is small enough), the viscosity should 
decrease with increase in M p due to the tempera- 
ture effect. For ATcmb = 10000K (Fig.[SH), F CM b 
is approximately proportional to M p . Even for 
ATcmb = 1000 K, Fcmb increases with M p for 
low mass regime (M p < 1M®). Thus, Fcmb has 
a peak at some value of M p for a given value of 
ATcmb- Figures [5}d and c show that the critical 
planet mass (M PjC ) at which Fomb takes the max- 
imum value is 2M® for AT C mb = 2000K and 5M® 
for 5000K. We empirically found that 



M 



AT, 



CMB 



1000K 



Ma 



(18) 



Since the mantle viscosity depends on the acti- 
vation volume, V* (Eq. [15]), and the values of 
V* may have uncertainty at high pressure, we 
also performed calculations with a smaller value 
of V*. Figures H] show the results with V* = 
3 x 10~ 6 m 3 mol -1 . Due to the weakened pres- 
sure effect, M p c is increased by a factor of a few. 
In Figures [HI the viscosity is artificially increased 
(Eq. [15] by a factor of 6000) in order to compen- 
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sate the smaller value of V* and reproduce Earth's 
observed values. Note that the artificial increase 
does not affect M PjC . 

The critical planetary mass M P]C is approxi- 
mately derived by the M p -dependence of the man- 
tle viscosity at CMB. Since we empirically found 
that P C mb ~ (Mp/M^P^cMB and T CMB ~ 
5AT CMB , Eq. ([15)1 is reduced to 



E* 



77cmb(M p ,T CM b) oc exp 



( M p \ 



PmV* 



tiRTcmb 



(19) 

When the argument of exponential is larger than 
unity, the viscosity rapidly increases with M p to 
depress Fcmb, because F C mb oc vj~ 1 / 3 . For M p < 
M PiC , we found that Fcmb M p . When the ar- 
gument exceeds some critical value (C > 1), the 
viscosity enhancement eventually overwhelms the 
factors for the positive M p -dependence of -Fcmb- 
Thus, M Pi c is given by the value of M p with which 
the argument of exponential is ~ C, 



M, 



p,c 



5nCRAT CM B - E* M ffi 

V* -P©CMB 

BnCRATcuB M e 



V* 



P 



V* 



C ATqmb 

10000A" I 10 x 10- 6 m 3 mol- 1 



©CMB 
1 



M ffi , (20) 



which explains the dependences on ATcmb and V* 
that we found numerically (If we adopt C ~ 10, 
the numerical factor is also explained.) 

3.4. Lifetime of intrinsic magnetic fields 

The lifetime of the intrinsic magnetic fields 
is calculated for AT CM b = 1000, 2000, 5000 and 
10000 K with a fixed value of V* = 10 x 
10 m mol . The results arc summarized in 
Fig. [TU] It is clearly shown that the lifetime de- 
clines for M p > M P:C by the increase in the mantle 
viscosity due to the pressure effect that we dis- 
cussed in details in the previous subsection. The 
results for V* = 3 x 10 — 6 m 3 mol — 1 show a similar 
property. 

In Fig. [TU1 the lifetime weakly increases with 
M p for M p < M PtC . The dependence is explained 
as follows. The lifetime is approximately given by 
Tnfc ~ ^/(ScmbFcmb), where E th is thermal 
energy of the core and Sqmb is surface area of 



the core. According to our calculation, Scmb oc 

1/2 2/3 

M p ' rather than M p ' due to self-compression. 
Figures 7 and 13 show that Fcmb M p and F t h oc 
Mp for a fixed Tqmb ■ As we mentioned in section 
3.3, T C mb ~ 5AXcmb- Thus, for a fixed AT C mb, 
it is predicted that r lifc - F th /(5'cMBFcMB) oc 



2-1/2-1 
P 



1/2 

Mp , which is consistent with the 



numerical results in Fig. 1101 

When Mp > M PiC , the higher mantle viscosity 
due to the effect of higher pressure significantly de- 
presses heat transfer at the bottom of the mantle. 
The suppressed heat flux cannot maintain the vig- 
orous core convection. As a result, the magnetic 
field lifetime is rather shorter for M p > M PjC . 

Figure [10] also shows that the lifetime for 
M p < Mp. c does not depend on ATcmb at all. 
The initial temperature is high enough to over- 
come the pressure-dependency even for the case 
of ATcmb = 1000 K, resulting in the effective 
self- regulation of Fcmb- As a result, the lifetime 
does not depend on the initial value of ATcmb- 

3.5. Strength of magnetic fields 

The strength of magnetic fields is as impor- 
tant as their lifetime to discuss habitability of the 
planets. Here we evaluate the strength of the 
magnetic fields, u s ing th e scaling law derived by 
Christensen et al. ( 20091). For planets with suffi- 
ciently rapid spins, Christensen et al. ( 20091 ) de- 
rived the magnetic strength at the core surface as 



r n t;,, 1 / 2 o 1 /6 z? 1 / 



3 

conv ) 



(21) 



where p c is average density of the core and 
fio is permeability. If the magnetic moment 
is dipol e-dominant and t he dipolc moment is 
oc 1/r 3 ( Gaidos et al.ll2010l ). the strength of mag- 
netic dipole at the planetary surface is B sul { = 
B c {r c /r p f. 

With F conv = F C MB-F cond , we calculated B sm[ 
from our simulation results. Figure QT] shows the 
calculated i? sur f at t = 5Gyr for various M p and 
AT C mb- The strength monotonically increases if 
the initial ATcmb is sufficiently high. The rela- 

1/3 

tively weak dependence (B oc M p ) comes from 

1/3 

the adopted the scaling law, B oc F conv , and the 
numerically obtained relation, Fcmb cx M p . If 
ATcmb is not high enough, the pressure effect 
is dominant and the strength is significantly sup- 
pressed for Mp > M p c . 
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4. Conclusion and Discussion 

We have developed a numerical model to sim- 
ulate thermal evolution of various-mass terrestrial 
planets in habitable zones. The density distri- 
bution of the planetary interior is calculated by 
Vinet EOS taking into account pressure depen- 
dence. Using the interior structure model, we cal- 
culate heat transfer through mantle, using the as- 
trophysical mixing length theory modified to man- 
tle convection. The modified mixing length theory 
is easily applied to multi-layer convection that may 
be dominated convection mode in super-Earths. 
We have calibrated the modified mixing length 
theory with the conventional parametrized con- 
vection model and the boundary layer theory, in 
simple one-layer convection cases. 

With nominal parameters of surface temper- 
ature T sur f = 300K, a mantle-core mass ration 
Cm/c = 7:3, initial core impurity xqs of 10 wt%, 
and initial temperature gap at CMB ATcmb = 
1000K, our model for M = 1M® reproduces sur- 
face heat flow and inner core radius of the present 
Earth. With different parameter values suitable 
for Mercury, Venus, and Mars, our model also re- 
produces the results that are not inconsistent with 
present magnetic activity of these planets. 

With this model, we calculated thermal evolu- 
tion of terrestrial planets with mass M p = 0.1- 
IOA/q in habitable zones, using the nominal pa- 
rameters, to study lifetime of intrinsic magnetic 
field that is one of the important factors for the 
planets to be habitable. We found from the nu- 
merical calculations that the lifetime is maximized 
at 



Mr, 



AT, 



CMB 



V* 



1000 A" \ 10 x lO-^mol" 1 



(22) 

where V* is activation volume of mantle material. 
Planets with smaller masses cool more rapidly, so 
that they cannot maintain core heat flux to gen- 
erate dynamo long enough. For M p > M PjC , the 
rapid increase in the mantle viscosity caused by 
high pressure significantly depresses heat trans- 
fer throughout the mantle and hence that in the 
core. As a result, dynamo cannot last long. Al- 
though the temperature effect tends to decrease 
the mantle viscosity as planetary mass becomes 
large, the pressure effect to increase the viscosity 
overwhelms the temperature effect for M p > M pc . 



With the numerically obtained empirical relation, 
Tqmb ~ 5 ATcmb, we can analytically derive 
Eq. (j2"2"|) from the Arrhenius-type formula for the 
mantle viscosity that we adopt (Eq. [15]). 

We found that while the lifetime of magnetic 
fields does not depend on ATcmb for M p < M PtC , 
it sensitively depends on ATcmb for M p > M p c 
because M PiC cx AT CM b (Eq. [22)- The initial 
ATcmb j that is, the initial temperature profile 
of planetary interior, is one of the most uncer- 
tain parameters, because it highly depends on 
the processes of planetary formation and differ- 
entiation of the planetary interior. As is shown 
by SPH simulations, if a planet undergoes gi- 
ant impacts, its metallic core is heated as high 
as several tens thousands K for M v <~ 1 M ffl 



( Canud 2004). On the other hand, if a planet 
accreted from small planetesimals without giant 
impacts, the initial temperature profile is deter- 
mined by the balance between gravitational en- 
ergy buried by planetesimals and thermal trans- 
fer efficiency through rocky mantle. The pro- 
cess includes crystallization of magma ocean and 
depends on the mechanical property of molten 



depends on tne m echanical proper ty oi molten 
mantle (lAbe fc Matsuil Il986l : IZahnle et al.lll988l 
Senshu et al.ll2002[) . Thus, to evaluate the lifetime 



of magnetic fields, in particular for super-Earths 
that arc likely satisfy M p > M PtC , detailed anal- 
ysis for accretion and early thermal evolution of 
terrestrial planets are needed. 

It is also found that higher initial tempera- 
ture profile delays the inner core nucleation. For 
super-Earths, in order to maintain magnetic field 
more than 10 Gyr, the initial temperature has 
to be high enough to overwhelm the pressure- 
dependence. However, in that case, the temper- 
ature of the core center never reaches its con- 
densation temperature and the inner core cannot 
grow. Some geo-dynamo simulations suggest that 
the presence of the inner core s tabilizes the dipole 
mom ent of geomagnetic field ( Sakuraba fc Konol 
1999). It is also suggested that because thermally 



driven convection is not sufficient to drive dynamo 
action against the ohmic dissipatio n within the 
core of Earth ( Gubbins et alJ [2003). the compo- 
sitional convection induced by light elements re- 
leased to the outer core by solidification of the in- 
ner c ore plays an essential role in dynamo genera - 
tion (jStevenson et al.lll983l : fGubbins et all 120041 ) ■ 
Since our results (Figures 6) show that inner 
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core is not nucleated and compositional convec- 
tion docs not occur for M p > 5M®, dipole mag- 
netic fields of super-Earths might not be stable. 

The existence of magnetic field of extrasolar 
planets could be directly detected by the polar- 
ization observation of the photon from transiting 
planets or detection of trapped by the mag- 
netic fields in the future. Another possibility of 
the detection of planetary magnetic field is, al- 
though it is indirect, observation of composition 
of planetary atmosphere or atmospheric tail. If 
the planet has intrinsic magnetic field, its atmo- 
sphere could keep H 2 molecules for long period. 
Venus may have lost H2O molecules o n a short 
time scale (IBullock fc Grinspoonl 120011) . Thus if 
water series molecules, such as H 2 0, H 3 0, and 
HO, were detected in the planetary atmosphere, 
it would indicate the existence of intrinsic mag- 
netic field, although super-Earths might be able 
to sustain the HxO molecules in the atmosphere 
by their high gravity even without the protection 
by magnetic fields. 

We need to elaborate our thermal evolution 
model, by considering details of mantle convection 
mode that is affected b y phase transition between 
7-spinel to perovskite (jChristensen fc Yuenlll985l ) 
at upper/lower mantle. We also should take into 
account further miner al transitions suggest ed by 
ab initio calculations ( Umemoto et al. 20061 ) that 
may appear in super-Earths, because they may 
affect internal density structure and the mantle 
convection mode. 

The abrupt enhancement in the mantle vis- 
cosity due to the pressure effect relies on the 
Arrhcnius-type formula for the mantle viscosity we 
adopt here. The critical mass beyond which the 
pressure effect dominates is inversely proportional 
to activation volume (Eq. [22]). Thus, detailed 
rhcological properties affect habitability of super- 
Earths. The values of the activation volume are 
not clear at such high pressure as in deep mantle in 
super-Earths. The mechanism to in hibit dynamo 
activit y in super-Earths proposed bv lGaidos et al 



also depends on high pressure material 
properties (melting and adiabatic curves), which 
also need to be confirmed. These provide new 
motivations to high pressure experiments and first 
principle simulations. Super-Earths provide good 
links between astronomy and high-pressure mate- 
rial science. 
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Appendix A. Comparison among Nu-Ra re- 
lationship model, thermal boundary layer 
model and mixing length theory model. 

In evaluation of thermal transfer of mantle con- 
vection, we_cj3m£are_Uie_mo^^ 
theory (jSasaki fe Nakazawalll986t lAbdll995l ) with 
the conventio nal parameterized conve ction model 



(PCM; e.g., ISharpe fc Peltierl Il979t ) and com- 



mon ly used thermal bound ary layer model (BLT; 
e.g.. lStevenson et"aL 1983 ). 

The original mixing leng th theory (MLT; e.g., 
Vitensd Il96l ISpiegell Il963r ) is often used in the 
thermal transfer within the st ellar interior to sim- 
ulate the stellar evolution. Sasaki fc Nakazawa 



()1986l ) modified the mixing length theory for very 
low Reynolds number convection in which the ver- 
tical flow is characterized by the Stokes velocity 
determined by a balance between buoyant force 
and resident force of viscosity rather than by free 
fall velocity. In the modified version, a distance 
(D) from the closest barrier such as the CMB or 
the top of the mantle layer is adopted for the mix- 
ing length £, while in the original theory, the den- 
sity scale height is usually adopted for £. 

The PCM uses the empirical Nu-Ra relation- 
ship, 



. Ra 

Nu = C[ 

1 Ra c 



1/3 



(23) 



where Nusselt number represents the ratio be- 
tween total heat flux and heat flux only due to 
conduction without convection, 



Nu 



-Ftotal 
-^cond 

AT 



(24) 



(25) 



and Raylcigh number is a dimcnsionlcss number 
representing the strength of convection, 



Ra = 



pgaATd 3 

KT] 



(26) 
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where g, a, AT and d are gravitational accelera- 
tion, thermal expansion, temperature difference 
between top and bottom and thickness of con- 
vective region, respectively, and Ra c is critical 
Rayleigh number (~ 650) for thermal convection. 
Because when Ra ~ Ra c , Nu must be ~ 1, £ 
is 0(1). Sotin et al. (1999) derived ( ~ 1.5-2.0 
through 3D fluid dynamical simulation although 
the value of £ is somewhat lower in high Ra re- 
gion. We here adopt £ = 1.7. 

From eqs. (|23|) to (|25|) . total heat flux through 
a fluid layer is represented by Rayleigh number as 

1/3 



F 



total 



— 

' \Ra c 



F r , 



(27) 



This model is very simple, but Ra is "mean" value 
of the whole mantle that is difficult to evaluate for 
real mantle in which viscosity changes by order 
of magnitude throughout the mantle. In particu- 
lar, it may not have enough resolution to evaluate 
Fcmb that we are concerned with in the present 
paper. 

In the BLT, heat flux is evaluated in the 
boundary layer. The thickness of boundary layer 
is estimated by an assumption that the layer 
is marginally stable against thermal instability. 
Then, the local Rayleigh number of the thermal 
boundary layer (Rai) is nearly equal to the critical 
Rayleigh number for thermal instability, that is, 



Rai = ( EEL) A 

\ K V / i 



TB 



Ra r 



(28) 



where S is thickness of thermal boundary layer and 
ATtb is temperature difference between the bot- 
tom and the top of the boundary layer, and sub- 
script "/" denotes the values in the thermal bound- 
ary layer. Thus, the heat flux through the layer is 
calculated as 

AT TB 



-Ftotal 



fc- 



wherc 6 is calculated as 



s = C 



O(l) 



(29) 



(30) 



is determined as follows. 

AT/2, so 



'(—) 


Ra c 


\pgaj 


i ATtb. 



and the factor (J 

If k, r\, p, g, and a are constant, ATtb ~ 
that 



(31) 



To be consistent with 3D fluid dynamical simula- 
tion by Sotin ct al. (1999), we set 2 4 / 3 C = 1/1.7, 
that is, £' = 0.23. 

Since the heat flux is expressed by quantities 
only in the thermal boundary layer (eqs. [27] and 
[28]), which is localized in the mantle, the BLT 
has better resolution than the PCM, in particular, 
for evaluation of -Fcmb- However, since the values 
of viscosity change by order of magnitude even in 
the thin thermal boundary layer, it is not clear 
which value has to be chosen as a representative 
value of the viscosity in Eq. [30] . For the terrestrial 
planets in our Solar system, observational data can 
be used to constrain the uncertainty 

Since the modified MLT uses local values of 
physical quantities (Eq. [13]), it quantitatively 
agrees with the BLT for wide range of parameters 
as shown below. There is no uncertainty for choice 
of a representative value of viscosity in the MLT, 
while choice of the mixing length has uncertainty. 
For calculation of thermal evolution of the Earth, 
it is proposed that a distance (D) from the closest 
barrier such as the CMB or the top of the man- 
tle layer is appropriate for £ j Sasaki & Nakazawa 
1986[ lAbeiri995l ISenshu et al.ll2002t iKimura et ~ 
20091 and references therein). Through compari- 
son with the calibrated PCM and BLT, we adopt 
£ = Q.82D as shown below. 

To compare these models, we calculate the heat 
flux in the case of radially constant r\ with the indi- 
vidual calibrated models. Internal heat generation 
due to radioactive elements is neglected. Figure 10 
shows the heat flux at the base of the mantle as 
a function of Ra, obtained by each model. The 
values are normalized by F con d, that is, Nusselt 
number. Although the MLT does not assume the 
relation of Nu cx Ra 1 / 3 , it produces the relation. 
To match the absolute values, we set £ = 0.82D. 
The maximum value of £ is proportional to d and 
heat flux is proportional to £ A (Eq. [14]) The sen- 
sitive dependence on d is canceled out to result in 
the rather weak dependence, Nu oc Ra 1 / 3 oc d, 
because we found that [(dT/dr) — (dT/dr) s ] de- 
creases with increase in £ (Eq. 11411. Analytical 
argument for it is found in I Abel (|l995l ). 

We also examined a case in which the viscosity 
is strongly temperature-dependent, 

r,(T) = m e X p[log(r/ 1 / % )(l - T)], (32) 
where t]q and r/i arc viscosity at the top (T = 0) 
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and the bottom (T = 1) of convective region. Fig- 
ure [13] shows Nusselt number obtained by PCM, 
BLT and MLT as a function of Ra. In the PCM, 
Ra is a mean value for a whole mantle. The repre- 
sentative viscosity is evaluated using average tem- 
perature of mantle, that is, T = 0.5 if mantle is 
thermally equilibrated because the PCM assume 
constant heat flux throughout mantle. In the BLT 
and MLT, the heat flux is evaluated by local quan- 
tities. The BLT and MLT produce the same heat 
flux within 1% in all cases, while the results by the 
PCM deviate from those by the BLT and the MLT 
for high Ra or high rjo/rji. These results show 
that MLT is as good as BLT to calculate ther- 
mal evolution of terrestrial planets. Since MLT 
is more easily to be applied for super-Earths that 
may have barriers for convection in their mantle 
(section 2.3), we adopt MLT. 

Appendix B. On the magnetism of planets 
in Solar system 

In order to confirm the validity of our model, 
we show that our model produces thermal evo- 
lution for individual terrestrial planets in the 
Solar system that is not inconsistent with their 
current magnetic activity, with appropriate non- 
nominal parame t er va lues, in a similar way to 



Stevenson et al.1 <| 19831 ) . Currently, Earth and 
Mercury have self-generating magnetic fields in- 
duced by dynamo action, while Venus and Mars 
do not (although some parts of the Martian 
crust have remnan t magnetic field in the past 



(jAcuna et al.lll999h ). 

To apply our model to Mercury, Venus and 
Mars, we need to use non- nominal parameter val- 
ues: 

• Mercury: Cm/e — 3 : 7 (a significantly 
large metallic core) and T sur f = 440K. 
These are observed values. We also tested 
sm aller values of % n s = .01,0.05 according 

We also tested 



to Stevenson et al 



(119830. 

higher mantle viscosity than Eq. (15) by 
multiplying viscosity increase factor A?; = 
100 

• Venus: T sur f = 737K, while the nominal val- 
ues are used for Cm/e and £os- We also tested 
higher mantle viscosity as well as in the case 
of Mercury. Note that two-layer convection 



model is used for Venus, because the spincl- 
perovskite transition also could work as a 
barrier for Venusian mantle. 

• Mars: T sur f = 210K. Cm/c is nominal value 
and xqs = 0.1,0.15,0.2. The standard for- 
mula, Eq. (15), is used for mantle viscosity. 

A?7 is viscosity increase factor due to lack of 
water in the case of Mercury and Venus. It is 
suggested by experiment that dry rock has factor 
of 100 higher viscosity than that of hydrated rocks. 
Thereby, we multiply Ary = 100 in the case of 
Venus and Mercury. 

The lifetime of magnetic fields calculated by our 
model is shown in Fig. 1141 In order to be consis- 
tent with current Mercury, Venus and Mars, the 
lifetime must be longer than 4.5Ga for Mercury 
and shorter than 4.5Ga for Venus and Mars. Be- 
cause Martian crust of age ~ 4Gyr retains pale- 
omagnetic field, the lifetime of Martian magnetic 
field may be longer than 0.5 Gyr. 

Figures [14] show that for Mercury, the lifetime 
is longer than 4.5Ga for relatively small values of 
.tqs (~ 0.01 — 0.05) except for extremely small 
ATcmb (< 200 - 300K). The relatively long life- 
time is resulted by nucleation of inner core due to 
lower solidification temperature corresponding to 
small values of xqs- If the nominal value of Xqs 
is used, the lifetime is short. The small value of 



c cnR fo r Mercury was discussed by iStevenson et al 
(jl983h . 



The predicted lifetime of magnetic field for 
Venus is quite short for relatively high mantle vis- 
cosity (Arj > 100). Observation suggests that 
Venus is lack of H 2 0. That may be due to run- 
away greenhouse effect of H 2 itself and conse- 
quent dissipation by UV dissociation and heating 
of the molecules. Because melting temperature of 
the mantle viscosity is lowered by H 2 0, relatively 
high mantle viscosity is more likely, although we 
do not know exact values of Venus' mantle viscos- 
ity. 

The predicted lifetime of magnetic field for for 
Mars is longer than 1 Gyr but shorter than 4.5Gyr, 
if initial AT CM b is ~ 10 - 500K. If Mars has 
never undergone giant impacts that cause signif- 
icant heating of metallic core, such low initial 
AT C mb is likely. 

Thus, with non-nominal parameters that reflect 
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distance from the Sun and accretion history of in- 
dividual planets, our model can produce the re- 
sults that are not inconsistent with the current 
terrestrial planets in the Solar system. However, 
in order to clarify intrinsic physics in generation 
of magnetic field of extra solar terrestrial planets, 
we focus on the results with the nominal parame- 
ters (T surf = 300K, Cm/c = 7:3, and x 0S = 0.1), 
which correspond to the parameters of terrestrial 
planets with the same compositions as the Earth 
in habitable zones. 

Appendix C. Energy in the core 

Figure [T5l shows the thermal and gravitational 
energy, released latent heat, and their sum as a 
function of Tomb for the nominal cases with plan- 
etary mass M = 1,2,5 and 10M®, which are cal- 
culated by the procedures in section 2.2. We set 
that each value is zero at the temperature at the 
initiation of inner core growth. As is shown in this 
figure, the loss of thermal energy occupies about 
one third of the total energy loss of the core for 
the case of M = 1M®. Released latent heat corre- 
sponds to about one fifth of the total energy loss, 
which depends on Tcmb because of the nonlinear 
density-dependency of the melting temperature of 
metal. 

The gradient of the total energy in Fig. [15] cor- 
responds to an effective specific heat of the core. 
The total heat capacity is twice larger than the 
specific heat of thermal energy solely just after the 
inner core initiation (Tcmb — 4100K), while their 
values converge as temperature decreases. This is 
because impurity concentration increases with the 
temperature decrease in the outer core. Inner core 
growth is moderated by the depression of melting 
temperature of outer core due to the concentra- 
tion of impurities into outer core. Released grav- 
itational energy and latent heat become smaller 
than thermal energy as the temperature decreases. 
Note that the gravitational energy released by the 
thermal contraction of the core also works as re- 
sistance to cooling of the core (see section [2~2]) . 

The ratio of gravitational energy, latent heat 
and thermal energy is varied with planetary mass. 
Thermal energy is more dominant than other en- 
ergies for more massive planet. It means that the 
gravitational energy and latent heat are not main 
energy source to drive dynamo action within cores 



of massive super-Earths. This is mainly because of 
the change in slope of adiabatic curve within core. 
The higher gravity causes steeper adiabatic ther- 
mal structure, and then core posses large amount 
of thermal energy inside it for the case of massive 
planets. This is also the reason why the effective 
specific heat of core is increased as planetary mass 
increases. 
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Table 1: Physical properties of mantle and core components we adopted (j Valencia et al.ll2007D 



material 


Po 


K 


K 


7o 


q 


0o 


Refs. 




(kgm" 3 ) 


(GPa) 












ol 


3347 


126.8 


4.274 


0.99 


2.1 


809 


a 


wd+rw 


3644 


174.5 


4.274 


1.20 


2.0 


908 


a 


pv+fmw 


4152 


223.6 


4.274 


1.48 


1.4 


1070 


a 


ppv+fmw 


4270 


233.6 


4.524 


1.68 


2.2 


1100 


b 


Fe 


8300 


164.8 


5.33 


1.36 


0.91 


998 


c,d 


FeS 


5330 


126 


4.8 


1.36 


0.91 


998 


c,d 



(Stixrude fe Lithgow-Bertellonil |2005h . b (|Tsuchiva et al.l |2004) . c (|Williams fe Knittld Il997h . d 



(jUchida et alj|2001l ) 



Table 2: Parameters of radiogenic elements we adopted (Ivan Schmudfl995h 



element 


(7( PP b) 


iiX/iWkg- 1 ) 


A(yr- i ; 




K 40 


28.0 


29.17 


5.54x10" 


10 


Th 232 


76.4 


26.38 


4.95x10" 


n 


TJ235 


0.14 


568.7 


9.85x10- 


10 


TJ238 


20.1 


94.65 


1.551x10 


-10 



Table 3: Parameter of the viscosity in upper mantle and lower mantle ( RanallilboOll ) 





S(Pa-"s- 1 ) 
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E*(10 3 Jmol- 1 ) 


V*(10- f5 m ;i mol- :L ) 


c(s- 


- 1 ) 


upper mantle 
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3.0 
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10 


10" 


-15 


lower mantle 


7.4xl0~ 17 


3.5 
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10 


10" 


-15 



Table 4: Physical property of upper mantle, lower mantle, and core (|Yukutakdl200m 





fc c (W mK- 1 ) 


C P (J kg-^K- 1 ) 




upper mantle 
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lower mantle 
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2.4xl0~ 5 


outer core 
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Fig. 1. — Radial density profiles for 0.1, 0.2, 0.5, 1, 2, 5, 10 M® planets (with inner core of 6 wt% of each 
core) obtained by our model in the nominal case. 
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Radius 



Fig. 2. — The schematic diagram of the procedure to obtain initial temperature distributions. For more 
detail, see text. 
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Fig. 3. — Time evolution of (a) the temperature profile, (b) the surface heat flux, (c) the inner core radius 
and (d) the heat flux through CMB for an Earth-mass planet with ATcmb = 1000-ftT. (a)Time evolution of 
internal temperature distribution in the case of ATcmb = 1000-K" with 1 M®. The planet radius of 6385km 
and core radius of 3375km. The surface heat flux declines to ~ 0.08 Wm" 2 and the inner core grows up to 
1200 km at 4.5 Gyr after, which is nearly equal to the present observed value of the inner core radius of the 
Earth. The core heat flux monotonically decreases with time but its time derivative discontinuously changes 
at the initiation of the inner core at around 2.2 Gyr. The solid curve represents the threshold flux (-Fcond) 
to maintain dynamo activity in the outer core (Eq. |17|). 
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Fig. 4. — Same figures as Figs. [3] except for two-layered mantle convection model with ATcmb =1000K. 
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Fig. 5. — Same figures as Figs. |3] except that two-layered mantle convection model with ATcmb — 800^ is 
considered instead of one-layered model with ATcmb =1000K. 
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Fig. 6. — Same figures as Figs. [3] except for M p = 5M®. 
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Fig. 7. Evolution of the core heat flux (left column) and inner core radius (right column) for M — (a)l, (b) 2, (c) 5, (d) and 10 

with various initial ATcmb ■ in th G left column, -F cr it for each ATcmb is expressed by thinner line with the same type. Some lines 
with different initial ATcmb arc overlapped by each other. For these parameters, the initial conditions do not affect the evolution since 
self-regulation of mantle heat transfer works due to temperature dependence of the mantle viscosity. Inner cores never nucleate in the 
cases of AT C mb = 2000, 5000, and 10000K for M p = 5M 9 and AT C mb = 5000 and 10000K for M p = 10A%. 
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Fig. 8.— Evolution of the core heat flux for AT CM b = (a)1000K, (b) 2000K, (c) 5000K, (d) 10000K with 
M p = 1, 2, 5, 10M®, respectively. 
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Fig. 10. — Lifetime of magnetic fields as a function of planetary mass (M p ) with various ATcmb for (a) 
V* = 10 X 10- 6 m 3 mol- 1 and (b) V* = 3 x Kr 6 m 3 mol _1 . 
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Fig. 11. — Strength of magnetic fields after 5 Gyr as a function of planetary mass (M p ) with various ATcmb 
for (a) V* = 10 x lO^i^mo!" 1 and (b) V* = 3 x 10- 6 m. 3 mol _1 . 



10° 



1CT 



10 u 



1CT 



10 u 



PCM n 

BLT 

MLT A 



B 



10 



10° 



Ra 



Fig. 12. — Temporal averaged Nusselt numbers obtained by PCM, BLT and MLT models. Open triangles, 
circles and squares represent the results of PCM, BLT and MLT models, respectively. 
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Fig. 13. — Temporal averaged Nussclt numbers obtained by PCM, BLT and MLT models. The viscosity is 
changed from the bottom to the top with the ranges of (a) Ar/= 10, (b) 10 2 , (c) 10 3 and (d) 10 4 , respectively 
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Fig. 14. — Lifetime of magnetic fields of (a) Mercury, (b) Venus and (c) Mars obtained through our simula- 
tions. The free parameters are initial temperature gap at CMB ATcmb and initial impurity concentration 
in core xo, viscosity increase factor A77, respectively. (a)£ m / c = 3:7, T sur f = 440K, and squares, circles, and 
triangles are corresponding to models with xo s = 0.01, 0.05, and 0.1, respectively. The filled and open symbols 
represent models with the standard viscosity and a higher viscosity multiplied by A77 = 100. (b)£ m / c = 7:3, 
Tsurf = 737K and x = 0.1. The filled and open symbols represent the same meaning as is in the case of 
Mercury. The square and triangles arc models in which one- and two-layered mantle convection are assumed. 
( c )Cm/c = 7:3, T surf = 210K. squares, triangles and circles are corresponding to models with x = 0.1, 0.15 
and 0.2. 
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Fig. 15. — Individual core energies as a function of X"cmb for the nominal case with M = 1, 2, 5, IOMq and 
•^os = 0.1. Dashed, dotted, and dot-dashed curves represent gravitational energy, latent heat, and thermal 
energy, respectively. Solid curve shows the total energy. The gradient of total energy corresponds to the 
specific heat of the core. 
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